PP01

Paquetes

library(spatstat)

## Loading required package: spatstat.data

## Loading required package: spatstat.geom

## spatstat.geom 3.2-9

## Loading required package: spatstat.random

## spatstat.random 3.2-3

## Loading required package: spatstat.explore

## Loading required package: nlme

## spatstat.explore 3.2-7

## Loading required package: spatstat.model

## Loading required package: rpart

## spatstat.model 3.2-11

## Loading required package: spatstat.linnet

## spatstat.linnet 3.1-5

## 
## spatstat 3.0-8 
## For an introduction to spatstat, type 'beginner'

Intensidad

Proceso Homogéneo

Intensidad promedio estimada se puede calcular de dos formas diferentes

npoints(swedishpines)/area(Window(swedishpines))

## [1] 0.007395833

intensity.ppp(swedishpines)

## [1] 0.007395833

intensity(swedishpines)

## [1] 0.007395833

Más detalles sobre el proceso puntual se pueden obtener al realizar

summary(swedishpines)

## Planar point pattern:  71 points
## Average intensity 0.007395833 points per square unit (one unit = 0.1 metres)
## 
## Coordinates are integers
## i.e. rounded to the nearest unit (one unit = 0.1 metres)
## 
## Window: rectangle = [0, 96] x [0, 100] units
## Window area = 9600 square units
## Unit of length: 0.1 metres

Si existen marcas, se calculará la intensidad para cada tipo

plot(amacrine)

intensity.ppp(amacrine)

##      off       on 
## 88.68302 94.92830

intensity(amacrine)

##      off       on 
## 88.68302 94.92830

Estimación intensidad

quadratcount(bei)

##            x
## y           [0,200) [200,400) [400,600) [600,800) [800,1e+03]
##   [400,500]     271       401       100        32         215
##   [300,400)     180       217       104        34         149
##   [200,300)     334         4        26        33         172
##   [100,200)     172        43        26       144          83
##   [0,100)       146        89       234       338          57

bei |> quadratcount( nx = 10, ny=10) |> plot(main = "Quadratcount")

Si queremos calcular la intensidad de acuerdo al conteo obtenido por el quadratcount (puntos dividido para el área), se puede usar la función intensity sobre un cuadrante.

bei |> quadratcount( nx = 10, ny=10) |> intensity(, image=TRUE) |> 
  plot(main = "Intensidad respecto a quadratcount")

Test de homogeneidad

Suponiendo una distribución de Poisson, la hipótesis nula asume que el proceso es un CRS.

qt <- unmark(amacrine) |> quadrat.test(7,7)
qt

## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  unmark(amacrine)
## X2 = 9.6667, df = 48, p-value = 8.386e-10
## alternative hypothesis: two.sided
## 
## Quadrats: 7 by 7 grid of tiles

plot(unmark(amacrine))
plot(qt, add = TRUE)

Comparar con la \(K\) función

env <- envelope(unmark(amacrine), Kest, nsim=19, rank=1, global=TRUE)

## Generating 19 simulations of CSR  ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 
## 19.
## 
## Done.

env

## Simultaneous critical envelopes for K(r)
## and observed value for 'unmark(amacrine)'
## Edge correction: "iso"
## Obtained from 19 simulations of CSR
## Envelope based on maximum deviation of K(r) from null value for CSR (known 
## exactly)
## Alternative: two.sided
## Significance level of simultaneous Monte Carlo test: 1/20 = 0.05
## ...........................................................
##      Math.label     Description                            
## r    r              distance argument r                    
## obs  hat(K)[obs](r) observed value of K(r) for data pattern
## theo K[theo](r)     theoretical value of K(r) for CSR      
## lo   hat(K)[lo](r)  lower critical boundary for K(r)       
## hi   hat(K)[hi](r)  upper critical boundary for K(r)       
## ...........................................................
## Default plot formula:  .~r
## where "." stands for 'obs', 'theo', 'hi', 'lo'
## Columns 'lo' and 'hi' will be plotted as shading (by default)
## Recommended range of argument r: [0, 0.25]
## Available range of argument r: [0, 0.25]
## Unit of length: 662 microns

plot(env)

bei |> 
envelope(Kest, nsim=19, rank=1, global=TRUE) |> 
  plot(main = "K para base bei")

## Generating 19 simulations of CSR  ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 
## 19.
## 
## Done.

Proceso no homogéneo

Se asume que la función de intensidad \(\lambda(u)\) depende de la ubicación \(s\) en el espacio.

Estimación no parámetrica

Los estimadores para la intensidad son:

• No correlacionada

\[ \tilde{\lambda}^{(0)}(u) = \sum_{i=1}^n \kappa(u-x_i)\]

• Correlacionada uniformemente \[ \tilde{\lambda}^{(U)}(u) = \dfrac{1}{e(u)}\sum_{i=1}^n \kappa(u-x_i)\]

• Corrección de Diggle \[ \tilde{\lambda}^{(U)}(u) = \sum_{i=1}^n \dfrac{1}{e(x_i)} \kappa(u-x_i)\]

para cualquier ubicación espacial \(u\) dentro de la ventana \(W\), donde \[ e(u) = \int_W \kappa(u-v) \]

y \(\kappa(u)\) es una función de densidad, es decir, \[ \kappa(u) \leq 0,\quad \forall u \in W \\ \int_{\mathbb{R}^2} \kappa(u) du = 1 \]

Observación

• Una selección usual para el kernel es considerar una distribución normal.

• La desviación estándar del kernel es el parámetro de suavizamiento o smoothing bandwidth. Un mayor valor de banda (bandwidth) incrementa el sesgo pero disminuye la varianza.

Teorema (Formula de Campbell)

Suponga \(f\) una función de ubicaciones espaciales \(u\), y considere la suma aleatoria \[ T = \sum_{i} f(x_i) \]

de los puntos \(x_i\) del proceso puntual \(X\). La formula de Cambell indica que \[ \mathbb{E}(T) = \mathbb{E}\left( \sum_{i} f(x_i) \right) = \int_{\mathbb{R}^2} f(u) \lambda(u) du \] donde \(\lambda\) es la función de intensidad de \(X\).

!!! pag 169, en la que se compara el uso de los estimadores.

Por defecto, se calcula el estimador uniforme del kernel, pero se puede indicar el estimador con la corrección de Diggle

amacrine |> density.ppp(diggle = TRUE) |> 
  plot(main = "Estimaición de la intensidad")

amacrine |> density.ppp(diggle = TRUE) |> 
  contour(main = "Estimaición de la intensidad")

amacrine |> density.ppp(diggle = TRUE) |> 
  persp(main = "Estimaición de la intensidad")

Selección de la banda de selección

El valor de la banda indica el nivel de suavizamiento que tendrá la intensidad. Si no se especifica la banda, se toma por defecto un octavo de la longitud del lado más corto de la ventana considerada. Sin embargo, se puede seleccionar el valor optimo considerando por ejemplo

b <- bw.diggle(amacrine)

## Warning: Berman-Diggle Cross-Validation criterion was minimised at right-hand
## end of interval [0, 0.0624]; use argument 'hmax' to specify a wider interval
## for bandwidth 'sigma'

b

##      sigma 
## 0.06237769

plot(b)

amacrine |> density.ppp(diggle = TRUE, sigma = b) |> 
  plot(main = "Estimaición de la intensidad con valor óptimo")

La banda puede variar bastante dependiendo del método que se utilice, donde cada uno tiene supuesto que pueden no adaptarse al modelo estudiado.

• bw.ppl asume un proceso de Poisson inhomogéneo.

• bw.diggle asume un proceso de Cox, que forma más grupos o clustes clustered que un proceso de Poisson.

• bw.scott

Métodos adaptativos

En todos los casos, hemos supuesto que tanto como la función de suavizamiento como la banda es la misma para todo el proceso. Eso puede mejorarse si se consideran métodos adaptativos. Referencia el paquete sparr en R.

Aquí f representa la fraacción de números del patrón puntual tomados de manera aleatoria y usados para realizar la partición de Dirichlet.

amacrine |> adaptive.density(f=1) |> plot()

amacrine |> adaptive.density(f=0.1, nrep = 30) |> plot()

## Computing 30 intensity estimates...1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
## 21, 22, 23, 24, 25, 26, 27, 28, 29, 
## 30.
## Done.

Otra opción es utilizar el vecino más cercano

amacrine |> nndensity(k = 10) |> plot()

Estimación paramétrica

• Intensidad homogénea: \(\lambda(u) = \lambda\)

• Intensidad homogénea por regiones \(\lambda(u) = b_j\) para todo \(u \in B_j\).

• Intensidad proporcional a un valor base: \(\lambda(u) = \theta b(u)\)

• Modelo exponencial: \(\lambda(u) = \kappa e^{\beta Z(u)} = \exp(\alpha + \beta(Z(u))\) donde \(Z(u)\) es una covariable espacial y \(\kappa\) y \(\beta\) son parámetros.

• Modelo de incidencia exponencial: \(\lambda(u) = b(u) \exp(\alpha + \beta Z(u))\)

• Modelo log lineal: \(\lambda_{\theta}(u) = \exp(B(u) + \boldsymbol{\theta}^\intercal \boldsymbol{Z}(u) ) = \exp(B(u) + \theta_1 Z_1(u) + \theta_2 Z_2(u) + \ldots + \theta_p Z_p(u) )\)

Ejemplo 1

Ubicación de árboles en una zona. Para el modelo log lineal

\[\lambda(u) = \exp(\beta_0 + \beta_1 S(u))\] donde \(S(u)\) es la pendiente en la ubicación \(u\).

fit <- ppm(bei ~ grad, data=bei.extra)
fit

## Nonstationary Poisson process
## Fitted to point pattern dataset 'bei'
## 
## Log intensity:  ~grad
## 
## Fitted trend coefficients:
## (Intercept)        grad 
##   -5.391053    5.026710 
## 
##              Estimate       S.E.   CI95.lo   CI95.hi Ztest      Zval
## (Intercept) -5.391053 0.03001787 -5.449887 -5.332219   *** -179.5948
## grad         5.026710 0.24534296  4.545847  5.507573   ***   20.4885

de donde \(\beta_0 \approx -5.39\) y \(\beta_1 \approx 5.03\).

Podemos observar el efecto que tiene la covariable en la intensidad

plot(effectfun(fit, "grad", se.fit=TRUE))

Otro modelo a considerar sería

\[\lambda(u) = \exp(\mu + \alpha_{E(u)} + \alpha_{G(u)} + \gamma_{E(u),G(u)})\] es decir,

fit <- ppm(bei ~ elev * grad, data=bei.extra)
fit

## Nonstationary Poisson process
## Fitted to point pattern dataset 'bei'
## 
## Log intensity:  ~elev * grad
## 
## Fitted trend coefficients:
##   (Intercept)          elev          grad     elev:grad 
##  -4.403426761  -0.007024363 -36.500458701   0.292813497 
## 
##                  Estimate        S.E.      CI95.lo       CI95.hi Ztest
## (Intercept)  -4.403426761 0.614763210  -5.60834051  -3.198513010   ***
## elev         -0.007024363 0.004192219  -0.01524096   0.001192235      
## grad        -36.500458701 5.261456400 -46.81272375 -26.188193651   ***
## elev:grad     0.292813497 0.036257125   0.22175084   0.363876155   ***
##                  Zval
## (Intercept) -7.162801
## elev        -1.675572
## grad        -6.937330
## elev:grad    8.076026

Podemos graficar también la estimación obtenida

plot(fit, how="image", se=FALSE)

Confidence intervals for the parameters

confint(fit, level=0.95)

##                    2.5 %        97.5 %
## (Intercept)  -5.60834051  -3.198513010
## elev         -0.01524096   0.001192235
## grad        -46.81272375 -26.188193651
## elev:grad     0.22175084   0.363876155

Ejemplo 2

Depositos de oro en Murchison, Australia. Primero pasamos los datos a kilómetros y calculamos la distancia a las fallas para cada punto.

mur <- solapply(murchison, rescale, s=1000, unitname="km")
dfault <- with(mur,distfun(faults))

estimados el modelo

\[\lambda(u) = \exp(\beta_0 + \beta_1 d(u))\] donde \(d(u)\) es la distancia más corta desde \(u\) a la falla geológica más cercana.

fit <- ppm(gold ~ dfault, data=mur)
fit

## Nonstationary Poisson process
## Fitted to point pattern dataset 'gold'
## 
## Log intensity:  ~dfault
## 
## Fitted trend coefficients:
## (Intercept)      dfault 
##  -4.3412775  -0.2607664 
## 
##               Estimate       S.E.    CI95.lo    CI95.hi Ztest      Zval
## (Intercept) -4.3412775 0.08556260 -4.5089771 -4.1735779   *** -50.73802
## dfault      -0.2607664 0.02018789 -0.3003339 -0.2211988   *** -12.91697

de donde \(\beta_0 \approx -4.34\) y \(\beta_1 \approx -0.26\).

plot(effectfun(fit, "dfault", se.fit=TRUE))

Comparar significancia de parámetros

Se compara con un modelos CRS y en este caso se rechaza.

fit0 <- ppm(bei ~ 1)
fit1 <- ppm(bei ~ grad, data=bei.extra)
anova(fit0, fit1, test="LR")

## Analysis of Deviance Table
## 
## Model 1: ~1   Poisson
## Model 2: ~grad    Poisson
##   Npar Df Deviance  Pr(>Chi)    
## 1    1                          
## 2    2  1   383.11 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Se estudia la inclusión o no de ciertas variables

fit2e <- ppm(bei ~ polynom(elev, 2), data=bei.extra)
fit2e1g <- update(fit2e, . ~ . + grad)
anova(fit2e, fit2e1g, test="LR")

## Analysis of Deviance Table
## 
## Model 1: ~elev + I(elev^2)    Poisson
## Model 2: ~elev + I(elev^2) + grad     Poisson
##   Npar Df Deviance  Pr(>Chi)    
## 1    3                          
## 2    4  1   419.24 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Selección de modelo mediante AIC

Preferimos el modelo con el AIC más pequeño.

fitprop <- ppm(bei ~ offset(log(grad)),data=bei.extra) 
fitnull <- ppm(bei ~1)
AIC(fitprop)

## [1] 42494.42

AIC(fitnull)

## [1] 42763.92

Comparación de residuos

diagnose.ppm(fit2e)

## Model diagnostics (raw residuals)
## Diagnostics available:
##  four-panel plot
##  mark plot 
##  smoothed residual field
##  x cumulative residuals
##  y cumulative residuals
##  sum of all residuals
## sum of raw residuals in entire window = 1.62e-09
## area of entire window = 5e+05
## quadrature area = 5e+05
## range of smoothed field =  [-0.005641, 0.006091]

Predicciones

fit <- ppm(bei ~ polynom(grad, elev, 2), data=bei.extra) 
lamhat <- predict(fit)
lamB <- predict(fit, locations=bei)

predict(fit) |> plot()

Simular el modelo ajustado

X <- simulate(fitprop);plot(X[[1]])

Imple

Dependencia con covariables

# Pesos

vols <- with(marks(finpines),
(pi/12) * height * (diameter/100)^2)
Dvol <- density(finpines, weights=vols, sigma=bw.ppl)
plot(Dvol)